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Stress intensity factors for interfaeial cracks in Fiber Metal Laminates (FML) are 
computed by using the displacement ratio method recently developed by Sun and Qian 
(1997, Int. J. Solids. Struct. 34, 2595-2609). Various FML configurations with single and 
multiple delaminations subjected to different loading conditions are investigated. The 
displacement ratio method requires the total energy release rate, bimaterial parameters, and 
relative crack surface displacements as input. Details of generating the energy release rates, 
defining bimaterial parameters with anisotropic elasticity, and selecting proper crack 
surface locations for obtaining relative crack surface displacements are discussed in the 
paper. Even though the individual energy release rates are noneonvergent, mesh-size- 
independent stress intensity factors can be obtained. This study also finds that the selection 
of reference length can affect the magnitudes and the mode mixity angles of the stress 
intensity factors; thus, it is important to report the reference length used with the calculated 
stress intensity factors. 


I. Introduction 

T HE NASA Aircraft Aging and Durability Program is currently supporting research efforts in the fracture and 
fatigue of Fiber Metal Laminates (FMLs). FMLs, such as Glare which contains glass fiber layers (about 0.005 
in/layer) and thin aluminum layers (about 0.016 in/layer), and CentrAL which contains Glare sublaminates in the 
middle and thicker outer aluminum layers * 1 * (0.32 in. or 0.16 in.), are being investigated. Schematics of the Glare and 
the CentrAL material systems are shown in Fig. 1. Note that GLARE® is a registered trademark of Structural 
Laminates Company, a joint venture of Aluminum Corporation of America (ALCOA) and Akzo Nobel to produce 
Glare; and CentrAL is a patent held by Alcoa, GTM Advanced Structures, and Delft University of Technology. 
Glare has been found to possess excellent fatigue and damage tolerance properties. Due to the presence of the fiber 
layers, the fatigue crack growth rate for Glare is much less than that for monolithic aluminum. Unfortunately, the 
manufacturing cost, related to milling, pretreatment, and storage of thin aluminum sheets and the labor intensive lay- 
up process, limits the configurations of Glare to approximately 6/5 lay-ups (6 metal sheets and 5 fiber prepreg layers 
in between). 1 Consequently, Glare is mainly used for fuselage panel applications, and cannot be used for wing box 
skins which require large material thickness. New CentrAL hybrid concepts 1 , which bond thicker outer aluminum 
layers to Glare, are being developed by Alcoa, Inc. and its partners to reduce manufacturing cost and satisfy fatigue 
and damage tolerance requirements of the thick lower wing skin panel applications. Research in fracture and fatigue 
of CentrAL FMLs is being conducted at NASA Langley Research Center under a Space Act Agreement between 
NASA and ALCOA. 

Delamination failure can significantly affect the fatigue and damage tolerance behaviors of FMLs. 
Delaminations often occur at the interfaces between metal layers and fiber layers. The crack growth behavior in an 
aluminum layer greatly depends on the delamination shape around (surrounding) the crack. 1 "’ To accurately predict 
the interfacial delamination growth in FMLs, an appropriate failure criterion needs to be established. 3 Due to the 
oscillatory characteristic of stress and displacement fields near the crack tip, the energy release rates of individual 
modes are not convergent with mesh refinement 4,5 thus, establishing the delamination failure criterion for FML 
interfacial cracks based on the individual mode of energy release rates may be problematic. On the other hand, the 
total energy release rate G T and the stress intensity factors of individual modes are well defined and are expected to 
be independent of mesh size. 4,h " s Thus, the ratios between the stress intensity factors, K u / K f and K ln / K t , and the 
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mode mixity angles, y/ n = tan '(K n / K , ) and i// l: , = tan ( K m / K ,) , are also expected to be independent of mesh 
size. Once the mode mixity angles are determined, the total energy release rate can be used to define the mixed 
mode failure criterion. 1,7 s 
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(a) Glare (b) CentrAL 

Fig. 1 Schematics of Glare and CentrAL cross-sections. 

The objective of this paper is to demonstrate that mesh-size-independent stress intensity factors are obtainable 
and can be efficiently computed for delamination cracks in FMLs. The stress intensity factors are computed with 
the displacement ratio method developed by Sun and Qian . 7 ’ 8 The inputs for the displacement ratio method are the 
total strain energy release rate, material parameters from anisotropic elasticity, and the near-tip relative crack surface 
displacements. 

This paper consists of four major sections. Section II describes background information related to the definition 
of the stress intensity factors, the explicit relationships between the near-tip stress fields and the stress intensity 
factors, and the explicit relationships between the displacement fields and the stress intensity factors. The Stroh 
formalism 9 of anisotropic elasticity used in establishing these relationships is given in the Appendix. The 
displacement ratio method for obtaining individual modes of stress intensity factors 7,8 is introduced next in this 
section. At the end of this section, a mixed mode fracture criterion is given to set the critical condition for an 
interfacial crack to propagate. In Section III, the stress intensity factors for delaminated FMLs subjected to the 
opening mode loading and the in-plane tensile loading are computed. Finite element models of FMLs with 
interfacial cracks are created and analyzed to provide the nodal forces and displacements needed for the Virtual 
Crack Closure Technique (VCCT ) 10 to obtain the energy release rates. The relative crack surface displacements and 
the total energy release rates are then used for the displacement ratio method to compute the stress intensity factors. 
Lessons learned such as how to avoid a potential pitfall in using the displacement ratio method and how to properly 
present the stress intensity factors are also given. Finally, Section IV gives concluding remarks and summarizes the 
findings of this study. 


II. Background Information 

In this section, relevant equations for computing the stress intensity factors for interfacial cracks in composite 
materials presented in Qian and Sun’s paper 8 are summarized. These equations are used for defining stress intensity 
factors, characterizing the near-tip stress and displacement fields, and presenting the displacement ratio method. All 
materials presented in this section can be applied for investigating the interfacial cracks in FMLs. 
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A. Stress Intensity Factors Definition 

The stress intensity factors definition for interfacial crack between two dissimilar anisotropic materials proposed 
by Hwu" and Wu 12 is used in this study. For generalized plane strain problems, such as the representative FML 
under an opening mode loading shown in Fig. 2, in which the strains are independent of the coordinate x, , the stress 
intensity factors can be expressed as 


K n 

K, 

K,„ 


■ = lim -JlirrA (^(r / r) IE " ^ A 1 • 


( 1 ) 


where the angular bracket represents a 3x3 diagonal matrix. s a (a =1,2,3) are constants related to the 
oscillation index s (see Appendix), and A is the eigenvector matrix associated with an eigenequation, Eq. (A 18), 
established by using the Stroh formalism 9 as discussed in the Appendix, f is an arbitrary chosen reference length 
parameter. 1 1-1 3 Stress intensity factors obtained for one reference length parameter can be easily converted to another 
reference length parameter 8 . For example, if the reference length parameter r is chosen to be the crack length 2 a , 
stress intensity factors for another r can be obtained by 

K(r) = A(((2 a/ry ie ° ))a^K(2 a) (2) 



Fig. 2 A representative fiber metaf faminate. 

The stress intensity factors defined in Eq. (1) have the following properties: 

1 . The units of the stress intensity factors, ( psi x -Jin ), are consistent with the conventional stress 
intensity factors. 

2. The stress intensity factors defined can be reduced to the classical stress intensity factors for a crack tip 
in homogeneous material. 

3. When the applied load is scaled by a factor, all intensity factors increase by the same factor. 

4. K n , Kj , and K m represent the intensity of singularity in the directions of x, , x 2 , and x 3 , 
respectively. 11 

5. Unlike the energy release rates computed by the VCCT, the stress intensity factors for interfacial 
cracks are convergent and not dependent on the virtual crack extension length (mesh size) A a . Due to 

the multiplication of ^(r/ r) w ° ^ in Eq. (1), the individual stress intensity factors are well defined 
and non-oscillatory, as evidenced by the explicit forms of the individual mode of stress intensity 
factors, which do not contain oscillatory terns similar to (Aa)^ (Refs. 4 and 8). 
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B. Characterizing Near Tip Fields 

Based on the stress intensity factors defined in Eq. (1), Hwu" found that the near-tip stresses and the near-tip 
relative crack surface displacements for an interfacial delamination crack can be expressed as 


and 
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Based on the material parameters presented in the Appendix, the near tip stresses in Eq. (3 ) for interfacial cracks in 
laminated composites can be explicitly written as 
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and the relative crack surface displacements in Eq. (4) can be explicitly written as 
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where W 2I and W 32 are elements of the W defined in Eq. (A17) in the Appendix. a f , b f , c f , and d f are also 
defined in the Appendix. g f , H x , H 2 , and complex constants c, , c, are defined as 

S f — d , W 32 , // — JV 2 X Kjj ~W 32 -Kjh , M 2 d f K u + ci j-K , 

«-('- ,; )“’*- (1 + 2 4‘osh(^) 
where s is the oscillation index defined in the Appendix. 

C. Displacement Ratio Method for Mode Separation 

The displacement ratio method 7,8 is based on the total energy release rate and the ratios of the near-tip crack 
surface displacements, obtained from finite element analysis results, to compute the stress intensity factors. The total 
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energy release rate can be computed by the VCCT. The following equation shows that G T can be expressed in 
terms of the stress intensity factors. 8 
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The total energy release rate equation is established by using the tip stresses and displacements in Eqs. (5) and (6) in 
Irwin’s crack closure integrals. 14 

The near tip relative crack surface displacement can be obtained from the same finite element analysis used for 
the VCCT. Based on Eq. (6), these crack surface displacements can be related to the stress intensity factors as 

Am, = c, K I + c 2 K h + c 2 K m 

Au 2 = d x Kj + d 2 K n + d 2 K m (8) 

Am 3 = e, Kj + e 2 K n + e 2 K m 

For brevity, the explicit expressions of coefficients c, , c/ ( , and e, (i = 1,2,3) are not presented here. From Eq. 8, the 
ratios K n / K 2 and K w / K ; can be determined by the displacement ratios Am, / Am, and Am 3 / Am, . Finally, these 
stress intensity factor ratios are used in Eq. 7 to obtain the K 2 , K u , and K m . The stress intensity factors can be 
obtained very efficiently by the displacement ratio method because the inputs for the displacement ratio method, the 
total energy release rate and the relative crack surface displacements, can be provided by a single finite element 
analysis. 

D. Failure Criterion 

The ratios of individual stress intensity factors, obtained by the displacement ratio method, are used to define the 
mode mixity angles y/ n and y/ n . Flere, tan y/ n =K u IK I and tan y/ u = K nl I K 2 . The critical condition for an 
interfacial crack to propagate can be given by a mixed mode fracture criterion 


A ( Vi.'. V ) A; (^12,1^13) 00 

where G c is the critical energy release rate (fracture toughness) that needs to be determined experimentally. Several 
specimens are available for this purpose such as the tensile sandwich test specimens 15 ' 1 ' 1 and the four point bend 
specimens. 17 y/ l2 and y/ u are the model mixity angles, defined using an arbitrary reference parameter r . Examples 
of y/ n and y/ u as a function of r are given in Section III.B. 

According to Ref. 18, the following steps may be taken to apply the mixed mode fracture criterion for assessing 
the interfacial crack failure in FMLs: (i) measure the fracture toughness G c of the interface crack as a function of 
mode mixity angle; (ii) calculate energy release rate G T and mode mixity for the interface crack in the FML 
component, and (iii) apply the fracture criterion to assess the failure load of the FML component. 


TIT. Stress Intensity Factors for Interfacial Cracks in FMLs 

FMLs with delamination cracks were analyzed to demonstrate that mesh-size-independent stress intensity factors 
can be easily and efficiently obtained by the displacement ratio method. Finite element analyses were performed to 
obtain the relative crack surface displacements and the nodal forces needed for using VCCT to compute the total 
energy release rate. These relative crack surface displacements and the total energy release rates are then used in the 
displacement ratio method to compute the stress intensity factors. Two FML configurations were used in this study. 
The first one, studied in Section III.A, is a representative fiber metal laminate with an interfacial delamination 
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between metal and fiber layers, and the second one, studied in Section III.B, is a CentrAL-like FML with multiple 
interfacial delaminations. 

The first FML configuration subjected to an opening mode loading was designed to test whether mesh size 
independent stress intensity factors can always be obtained for interfacial cracks in FMLs. Plane strain analyses 
using ABAQUS' 2D elements were performed for models with principal material directions aligned with the 
specimen’s x x and x 3 axes (see Fig. 3) and generalized plane strain analyses using ABAQUS' 3D elements were 
performed for models with principal material directions not aligned with the axes. 

The second FML configuration subjected to an in-plane tensile loading was designed to investigate the variations 
of stress intensity factors with delamination locations in a CentrAL-like FML with multiple interfacial delaminations 
under an in-plane tensile loading. The effects of the reference length parameter in stress intensity factors and mode 
mixity angles were also investigated. These computed stress intensity factors and mode mixity angles could be used 
in the interfacial failure criterion in Eq. (9) to determine the interfacial crack propagations. 

Material properties for the fiber and the aluminum layers of the FMLs are listed in Table 1. These properties are 
obtained from Ref. 19. 


Table 1 Material Properties of Unidirectional Glass/Epoxy and Aluminum Layers. 


Glass/Epoxy 

Aluminum 

E x (msi) 

7.977 

E (msi) 

10.5 

E 2 (msi) 

1.378 

V 

0.33 

E 3 (msi) 

1.378 



Ua 

0.33 



U 3 

0.33 



^23 

0.45 



G 12 (msi) 

0.8 



G, 3 (msi) 

0.8 



G 23 (msi) 

0.435 




A. Representative Fiber Metal Laminate under Opening Mode Loading 

A representative fiber metal laminate with a delamination under an opening mode loading shown in Fig. 2 was 
analyzed to demonstrate that mesh-size-independent stress intensity factors can be computed using the displacement 
ratio method. The delamination is located at the interface between the upper aluminum part and the lower four-layer 
Glass/Epoxy part. The lay-ups of the Glass/Epoxy part include [0 4 ], [+15/-15] s , [+30/-30] s , [+45/-45] s , [+60/-60] s , 
[+75/-75]s, and [90 4 ]. Only linear analyses were performed in this study. The thickness of the aluminum is 2.0 in. 
and the thickness of each fiber layer is 0.5 in. to assure high bending stiffnesses for both crack flanks to avoid 
geometrically nonlinear deformations. The delamination crack length is 10 in. and an equal and opposite vertical 
load P of 1 00 lb/in (distributed load in x, direction) is applied at the left end of each crack flank. The right end of the 
model is fixed. Plane strain and generalized plane strain finite element models, using ABAQUS® CPE8 and C3D8I 
elements, respectively, were created and analyzed with ABAQUS® /Standard. 20 Plane strain models using 2D CPE8 
elements were created for FMLs with [0 4 ] and [90 4 ] Glass/Epoxy lay-ups, and generalized plane strain models using 
3D C3D8I elements were created for all seven types of Glass/Epoxy layups to properly input the off-axis 
Glass/Epoxy material properties. The VCCT is used to obtain the energy release rates. The total energy release rate, 
obtained by summing the individual mode of energy release rates, is then used in the displacement ratio method for 
calculating the stress intensity factors. 
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1. Plane Strain Finite Element Analyses 

Four plane strain finite element models, representing three different crack tip mesh sizes, au = 0.0125, 0.025, 
0.05, and 0.1 in., were created (see the inset in Fig. 3). Two ABAQUS finite element analyses were performed. 
The first analysis assumed the glass fibers were laid up in a [O4] stack, and the second analysis assumed the glass 
fibers were laid up in a [90 4 ] stack. 


Aluminum 



Glass/Epoxy 


u 


V 


— ►! f* — A a 

Fig. 3 Typical plane strain finite element model. 


0 

0 


The finite element model with fiber in the 0-deg direction is shown in Fig. 3. The energy release rates, obtained 
by the VCCT, for the models with glass fibers in the 0-deg direction are plotted in Fig. 4. The opening mode energy 
release rate G ; and the sliding mode energy release rate G u are nonconvergent with increased mesh refinement, but 
the total energy release rate G T is well defined and independent of the crack tip mesh size (i.e. decreasing a a ). 

The individual stress intensity factors for the models with various a a were computed with the displacement 
ratio method, using two different reference length parameters (r =2 and 20 in). Note that the relative crack surface 
displacements for the displacement ratio method are taken at one element away from the crack tip. Figure 5 shows 
that the individual stress intensity factors computed are nearly constant for crack tip mesh sizes a a > 0.025 in. The 
stress intensity factor K u at ali = 0.0125 in. is less than those of larger mesh sizes. This may be due to the fact that 
the oscillatory displacement fields near the crack tip cannot be accurately modeled with the quadratic ABAQUS' 
CPE8 elements. It is reasonable to hypothesize that the relative crack surface displacements at a location further 
away from the crack tip may be better predicted by ABAQUS '. To test this hypothesis, the relative crack opening 
displacements were taken at the location 0.025 in. (two elements) away from the crack tip for computing the stress 
intensity factors for the model of ao= 0.0125 in. The stress intensity values of the smallest mesh size model are 
increased to be in line with those of the larger mesh size models as shown in Fig. 6. This indicates that mesh-size- 
independent stress intensity factors can be obtained if the relative crack opening displacements are not taken very 
near the crack tip. 

Furthermore, Figures 5 and 6 show that changing the reference length from f =2 to r= 20 in. can result in 
changes of both K , and Rvalues. For this case, the changes are non-proportional, a larger reduction occurs in the 
K u values than in the values. The non-proportional changes indicate that changing the reference length not only 
changes the magnitudes of K r and K u but also their mode mixity angle, t// p = tan ~'{K U / Kj) . Fortunately, the 
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Fig. 4 Energy release rates for AL/(0-deg Glass/Epoxy) laminate. 
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Fig. 5 Stress intensity factors for AL/(0-deg Glass/Epoxy) laminate with relative crack-tip 
displacements at one element away from the crack tip. 
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Fig. 6 Stress intensity factors for AL/(0-deg Glass/Epoxy) laminate (leftmost data computed 
with relative crack-tip displacements at two elements away from the crack tip). 


stress intensity factors computed by using one reference length can easily be converted to another reference length 
with Eq. (2). 

The energy release rates computed for the models with glass fibers in the 90-deg direction are plotted in Fig. 7. 
G r is well defined, however, G t and G u are nonconvergent with mesh refinement and can cross each other. This 
indicates that using G t and G n in an interfacial crack failure criterion may be problematic. The individual stress 
intensity factors for models with various aa were obtained by the displacement ratio method. Figure 8 shows that 
the individual stress intensity factors are independent of mesh sizes for aa > 0.025 in. The stress intensity factor 
K n at aa= 0.0125 in. is less than those of the larger mesh size models. As discussed for the 0-deg case, the 
oscillatory displacement fields very near the crack-tip may not be accurately predicted using the ABAQUS' CPE8 
elements. Again, by taking the relative crack surface displacements at a distance of 0.025 in. (two elements) away 
from the crack tip, the values of the stress intensity factors for the smallest mesh size model are nearly the same as 
those of the larger mesh size models as shown in Fig. 9. Note that the stress intensity factors for two reference 
length parameters are shown in Figs. 8 and 9. Again, Figs. 8 and 9 show that changing the reference length from 
f =2 to r =20 in. can result in the changes of the K , and K n values and the mode mixity angle, 

Yn =tan -\K II !K I ). 
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a a (in) 

Fig. 7 Energy release rates for AL/(90-deg Glass/Epoxy) laminate. 



a a (in) 

Fig. 8 Stress intensity factors for AL/(90-deg Glass/Epoxy) laminate with relative crack surface 
displacements at one element away from the crack tip. 
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Fig. 9 Stress intensity factors for AL/(90-deg Glass/Epoxy) laminate (leftmost data computed 
with relative crack-tip displacements at two elements away from the crack tip). 


2. Generalized Plane Strain Finite Element Analyses 

A generalized plane strain finite element model for the FML shown in Fig. 2 was created and analyzed for 
studying the effects of varying the fiber angles of the Glass/Epoxy lay-ups on the energy release rates and the stress 
intensity factors. A modeling technique that can enforce the generalized plane strain deformations in theXj-x 2 
plane is used. The finite element model, created by extruding the mesh pattern in Fig. 3 along the x 3 -direction, 
contains only one layer of brick elements, 0.05 in. in the x 3 -direction. The 3D ABAQUS® C3D81 finite element 
model has element size A a =0.025 in. around the crack tip. To ensure no variation of strains in the x 3 -direction, the 
corresponding nodes in the two parallel surfaces perpendicular to the x 3 -direction are constrained to have the same 
x 3 -displacements. The loading conditions and boundary conditions are the same as those shown in Fig. 3. 

Seven finite element analyses were performed for the models with different lay-ups of the Glass/Epoxy part. 
These lay-ups are [0 4 ], [+15/-15] s , [+30/-30] s , [+45/-45J s , [+60/-60] s , [+75/-75] s , and [90 4 ]. Because of the use of 
3D solid elements, the engineering constants associated with the three principal material directions of each 
orthotropic ply and the ply orientation can be directly assigned to each element. The energy release rates versus lay- 
up angles for the interfacial crack are plotted in Fig. 10. The values of G I , G n , and G T increase monotonically due 
to the increasing of flexibility of the Glass/Epoxy part. Note all G m values are near zero. The total energy release 
rates shown in Fig. 10 were used in the displacement ratio method to obtain the stress intensity factors. Figure 11 
plots the Kj, K u . and K IU versus lay-up angles, which represent the various lay-ups. For example, angle (+/-45 
Deg) on the abscissa is the [+45/-45J s lay-up. K , increases as the fiber angle increases, and the curve becomes flat 
for fiber angles greater than 60-deg. The maximum values of K n and K IU occur at the 45-deg angle. Note that for 
the same lay-up, a zero G m value in Fig. 10 does not correspond to a zero K m value in Fig. 11, indicating that there 
is not a one to one relationship between K U1 and G m . Each stress intensity factor is coupled with all three 
individual modes of strain energy release rates. K n becomes positive for the fiber angles greater than 25-deg and 
nearly constant for the fiber angles greater than 45-deg. The mode mixity angles can then be determined with the 
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computed K t , K n , and K nl values and be used in Eq. (9) for predicting the mixed mode delamination failure at the 
interface between the aluminum part and the Glass/Epoxy part shown in Fig. 2. 
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Fig. 10 Energy release rates for FMLs with four symmetric fiber layers. 
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Fig. 11 Stress intensity factors for FMFs with four symmetric fiber layers. 
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B. FML with Multiple Interfacial Cracks under In-Plane Tension Loading 

A CentrAL-like configuration with fractured aluminum layers and multiple delaminations under in-plane tension 
loading shown in Fig. 12a was analyzed. This configuration is selected to demonstrate that the analysis technique 
developed here can also be applied to predict the failure of multiple interfacial delaminations in the CentrAL FMLs. 
Note that all the fiber layers are intact. Figure 12b shows a cross-section of 0.1 16 in. in length and 0.176 in. in height 
in the x, - x 2 plane. The two outer aluminum layers are twice as thick as the internal aluminum layers. The 
thicknesses of both the aluminum and Glass/Epoxy layers are shown in the figure. Near the cracks in the aluminum 
layers, short delaminations exist between the aluminum and the Glass/Epoxy layers. This configuration was loaded 
with a uniform in-plane tension loading that can induce shear deformations at the interfacial crack tips. By 
exploiting symmetry, a finite element model of the upper right quadrant was created and analyzed. Figure 13 shows 
the symmetry boundary conditions, u=0 at the left ends of the Glass/Epoxy layers and v=0 at the bottom edge of the 
model. A uniform displacement, u=0.0015 in. was applied at the right edge. This upper right quadrant contains four 
delamination cracks and was modeled with ABAQUS 8-node plane strain, CPE8, elements. The strain energy 
release rates and the stress intensity factors were computed for each delamination crack. 

The deformed configuration shown in Fig. 13 reveals that the crack opening displacements are gradually 
decreasing from the outermost crack. Crack #1, to the innermost crack. Crack #4. This indicates that Crck #1 may 
have the highest Mode I energy release rate and stress intensity factor. The values of the energy release rates Gj, Gn 
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Fig. 12 Fiber metal laminate with multiple cracked aluminum layers and delaminations. 
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, and Gt for each crack are plotted in Fig. 14. The mode separation method presented in Section II. C is used to 
compute the stress intensity factors for each delamination crack. Two different reference length parameters, f = 
0.032 in. and 2.0 in., are used for computing the stress intensity factors to investigate whether the computed stress 
intensity factors and the mode mixity angles are dependent on the reference parameter used. The values of the stress 
intensity factors are plotted in Figs. 15a and 15b, and the values of the mode mixity angles are plotted in Fig. 15c. 
Figures 14 and 15 show that Mode-11 is the dominant mode under the in-plane tensile load, G n » G , and K n » K, . 
The Mode-I energy release rate data in Fig. 14 and the stress intensity factor data in Fig. 15a both show that the 
greatest values occur at Crack #1. 

The total energy release rates and mode mixity angles for the delamination cracks, obtained by using the two 
different reference length parameters, are shown in Fig. 16. These total energy release rates and mode mixity angles 
can be used with the mixed mode failure criterion, Eq. (9), to determine delamination failure. For the same crack, 
the mode mixity angle obtained by using one reference length parameter can be quite different from using another 
reference length parameter as shown in Figs. 15 and 16. As a result, it is necessary to provide the reference length 
parameter with the stress intensity factors reported. The user can convert the stress intensity factors to other 
reference length parameters by using Eq. (2). 
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Fig. 13 Deformed shape and boundary conditions. 



12 3 4 

Crack Number 


Fig. 14 Energy release rates for multiple cracks in the FML. 
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Fig. 15 Stress intensity factors and mode mixity angles for a multiple cracked FML. 
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Fig. 16 Total energy release rates and mode mixity angles for multiple cracks in FML. 
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Fig. 17 Stress intensity factors predicted by using relative crack surface displacements 
at various distances from the crack tip. 


16 

American Institute of Aeronautics and Astronautics 


Mode Mixity Angle (deg) 


The discussion in Section III. A. 1 notes that the stress intensity factors obtained by the displacement ratio method 
may be inaccurate if the relative crack surface displacements are taken at a location very near the crack tip. It can be 
useful to understand the sensitivity between the predicted stress intensity factors and the locations (distances from 
the crack tip) for taking the relative crack surface displacements. Relative crack surface displacements at four 
different locations taken from the finite element analysis results were used in the displacement methods to obtain the 
stress intensity factors. These locations are at 3.125%, 6.25%, 12.55% and 37.5% of the delamination crack length 
from the crack tip. Note that the delamination crack length in the finite element model is 0.016 in. (see Figs. 12). 
Stress intensity factors computed with these relative crack surface displacements for four different locations are 
plotted in Fig. 17. The stress intensity factors and mode mixity angles obtained by using the relative crack surface 
displacements at the first three locations are similar. Using the relative crack surface displacements at 37.5% of the 
delamination length, the stress intensity factors are all less than 3% from the previous values. This may indicate that 
the relative crack surface displacements can be taken at a distance of many elements away from the crack tip, and 
the computed stress intensity factors are insensitive to the location where the relative crack surface displacements 
are taken. 


IV. Concluding Remarks 

A procedure for obtaining stress intensity factors based on the total strain energy release rate, material 
parameters from anisotropic elasticity, and the near-tip relative crack surface displacement ratios was implemented 
for calculating the stress intensity factors of interfacial cracks in FMLs. This study showed that mesh-size- 
independent stress intensity factors are obtainable if the relative crack surface displacements are not taken at a 
location very near the crack tip. By using these stress intensity factors, mesh-size-independent mode mixity angles 
can be obtained and used with the experimentally determined critical total energy release rate for defining a mixed 
mode fracture criterion. The stress intensity factors can be obtained very efficiently by the displacement method 
because its inputs, the total energy release rate and the relative crack surface displacements, can be provided by a 
single finite element analysis. 

Various FML configurations with single or multiple delaminations and different loading conditions were 
investigated and presented in the paper. This study found that fiber orientations can significantly affect the Mode-1 
stress intensity factors for a FML configuration with an opening mode loading. For a CentrAL-like configuration 
with fractured aluminum layers and multiple delaminations, both the K k and K u stress intensity factors have their 
highest values at the outermost crack and the lowest values at the innermost crack, gradually decreasing the values 
of stress intensity factors from the outermost crack toward the innermost crack. 

This study also found that the selection of reference length can affect the magnitudes and the mode mixity angles 
of the stress intensity factors; thus, it is important to report the reference length used for calculating the reported 
stress intensity factors. The last finding is that the relative crack surface displacements can be taken at a distance of 
many elements away from the crack tip. The stress intensity factors computed by the displacement ratio method are 
insensitive to the locations where the relative crack surface displacements were taken. 

Appendix - Stroh Formalism, Bimaterial Eigenequation, and Bimaterial Parameters 

This appendix provides some background information related to how the eigenvalues and eigenvectors used in 
Eq. (1) are obtained. The eigenequation was derived by Hwu 21 based on the sextic formalism of Stroh. 9 A brief 
introduction of Stroh formalism is included here. For a bimaterial interfacial crack, Stroh formalism was used to 
establish the impedance matrices 21 M for the upper part material. Material #1, and the lower part material, Material 
#2, that define the crack plane. These impedance matrices were then used to form the bimaterial matrix M* that was 
used to establish the bimaterial eigenequation. 21 The forms of the matrices M and M will be given later. 

The stress-strain relations of an anisotropic system are given by 

Vij = C ijks U k,s ( A1 ) 


where cr. are the stresses, u k s are the strains, and C jjks are the stiffnesses. The equations of equilibrium are given 


as 


C ijks U k.s ij - 0 


(A2) 
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For a generalized plane strain deformation in which u k (k = 1,2,3) only depend on x, and x 2 of a fixed 
rectangular coordinate system, the displacements can be assumed as 9 ’ 22 ’ 23 

u k= a kf( z ) (A3) 

where z = x, + px 2 . p and a k are unknown constants to be determined. The equations of equilibrium can be 
rewritten as 


c i]ks + pSj2 ) (^1 + pS s2 )a k =0 (A4) 

or 

{Qi k\ + P(C nk2 +C i2 kl ) + p C ak2 }a k = 0 (A5) 


rl.O, i = j 

where S = 1 denotes the Kroneker delta. This can be written in matrix notation as 

v L 0.0, i * j 


Q + p(R + R r ) + p 2 T}« = 0 


(A6) 


where Q,R, and T are 3x3 matrices whose components are 

Q* = C im , R ik = c nk2 , T ik = C i2k2 {i,k = 1,2,3 ) ( A7) 

Eq. (A6) can be rewritten as 

(Q + pR)a = -p{ R r + pT)a (A8) 

By introducing a new vector b = (R r + pT)a and using Eq. (A8), the following expression can be derived. 


b = (R r + pT)a = 


-i(Q + ^R)a 

P 


Equation (A8) can be recast as a standard eigenvalue problem, 


N $ = A 


in which 


N = 


-T R r 
RTRAq 


T “1 

-RT 7 


, and £ = 



(A9) 


(A10) 


(All) 


Six eigenvalues of p and six associated eigenvectors of £ can be obtained. It can be shown that p cannot be real if 
the strain energy is positive. 2 ' Therefore, three pairs of complex conjugate p and three pairs of complex conjugate C 
are as follows 


P m = Pi and £ +3 =£. O' = 1,2,3) 


(A12) 
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Note that p t ( i = 1, 3 ) have positive imaginaiy parts and their associated eigenvectors, £ ii = 


, are used to define 


two matrices, A and B. 

A = [a, a 2 a 3 ],B = [b, b 2 b 3 ] 

The bimaterial matrix M* (Ref. 21) is defined with the impedance matrices, M. (i=l,2) as 

M* =M; 1 +M 2 ' 1 =i(A 1 B;‘-A 2 B 2 ) = -i(W + iD) (A13) 

where W and D are real matices and are related to Barnet-Lothe’s tensors S and L (Ref. 24) as 

W = S,L 1 -S 2 L 2 , D = L,' + L 2 (A14) 

In this study, the upper part, MAT#1, is isotropic, the A, B, 1 term may be computed directly from the explicit 
expressions of S[ and L 2 (Ref. 22) to avoid the problems related to the degenerate material by writing the A^j 1 


where 


AjB ," 1 = -(SjL ) 1 +iL)‘) 


Sj 


1 


l-2v 
2(1 -v) 


0 

1 

0 


-1 

0 

0 


> L 1 =M 


1 - 1 / 

0 

0 


0 0 


0 1 


(A 15) 


(Alb) 


v is the poisson’s ratio, and ju is the shear modulus. 

For FMLs, the D and W can be expressed as the following forms 8 



D n 

0 

0,3 ' 


" 0 

-W 

r '2l 

0 

D = 

0 

D 22 

0 

, w = 

W 21 

0 

-W 

vv 32 


. 0,3 

0 

033 _ 


0 

w 

” 32 

0 


(A17) 


The A in Eq. (1) is obtained from the following bimaterial eigenequation formulated by Hwu. 21 

(M* +e 2i ' r5 M*) A = 0 


(A18) 


The explicit solution of 8 was given by Ting 5 as 


8 a =-- + is a , a = 1,2,3 


(A19) 


where 


1 , l + fi n 

£•,=£• = In , £•,=-£, £•, = 0 , 

1 2;r 1 -P 


(A20) 
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and 


P = 



frfWD ' 1 ) 2 


(A21) 


where tr stands for the trace of a matrix and s is called the oscillation index. The three eigenvectors of Eq. (A 18) 
are A = [k, k 2 >.,] and are normalized by A T DA = I . Using the W and D forms shown in Eq. (A17), the /? and 
A can be explicitly expressed as 8 


ft = . 


\W 2 i cif + W n d f 
I D 22 e f 


A = 


-ia ( !b f c f ia f lb f c f W 32 / ^ Jb f j3D 2 . 

1 / c f 1 lc f 0 

id f / b f c f -id f / b f c f W 2I / ^b f j3D 2 . 


(A22) 


(A23) 


where the coefficients a f , b f , c f , d f , and e f are given as 


df — W 2i D 33 + W 32 D l 3 
b f = ft (D lt D 33 - 0,3 ) 

c f = t]2D 22 (A24) 

d f = W 32 D n +W 21 D I3 

e f=D n D 33 -D? 3 

The coefficients a f , b f , c . , d f , and e f , oscillation index s , and components of W and D matrices are designated 
as bimaterial parameters which are used in Eqs. (1) to (8). 
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